{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:58:06.468311",
     "start_time": "2016-08-17T07:58:04.883583"
    },
    "autoscroll": "json-false",
    "collapsed": false,
    "ein.tags": [
     "worksheet-0"
    ]
   },
   "outputs": [],
   "source": [
    "%load_ext autoreload\n",
    "%autoreload 2\n",
    "\n",
    "import sys \n",
    "from os import getcwd, path\n",
    "sys.path.append(path.dirname(getcwd()))\n",
    "import matplotlib.pyplot as plt\n",
    "import seaborn as sb\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "import pprint\n",
    "from cohorts.functions import *\n",
    "import lifelines as ll\n",
    "import patsy\n",
    "import functools\n",
    "import survivalstan\n",
    "from cohorts.utils import strip_column_names\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## set seeds for stan & rngs, to aid in reproducibility\n",
    "## (note: only reproducible within a machine, not across machines)\n",
    "\n",
    "seed = 91038753\n",
    "import random\n",
    "random.seed(seed)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:58:06.519332",
     "start_time": "2016-08-17T07:58:06.499244"
    },
    "autoscroll": "json-false",
    "collapsed": false,
    "ein.tags": [
     "worksheet-0"
    ]
   },
   "outputs": [],
   "source": [
    "from utils import data\n",
    "from utils import paper\n",
    "from utils.extra_functions import *"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# prep data "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## load data"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:58:11.840853",
     "start_time": "2016-08-17T07:58:06.520959"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inner join with tcr_peripheral_a: 25 to 25 rows\n",
      "inner join with ensembl_coverage: 25 to 25 rows\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inner join with tcr_peripheral_a: 25 to 25 rows\n",
      "inner join with ensembl_coverage: 25 to 25 rows\n",
      "{'dataframe_hash': 7698303973572390439,\n",
      " 'provenance_file_summary': {u'cohorts': u'0.4.0+3.gda968fb',\n",
      "                             u'isovar': u'0.0.6',\n",
      "                             u'mhctools': u'0.3.0',\n",
      "                             u'numpy': u'1.11.1',\n",
      "                             u'pandas': u'0.18.1',\n",
      "                             u'pyensembl': u'1.0.3',\n",
      "                             u'scipy': u'0.18.1',\n",
      "                             u'topiary': u'0.1.0',\n",
      "                             u'varcode': u'0.5.10'}}\n"
     ]
    }
   ],
   "source": [
    "cohort = data.init_cohort(join_with=[\"ensembl_coverage\",\"tcr_peripheral_a\"])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:21.487992",
     "start_time": "2016-08-17T07:58:11.931783"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inner join with tcr_peripheral_a: 25 to 25 rows\n",
      "inner join with ensembl_coverage: 25 to 25 rows\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "inner join with ensembl_coverage: 25 to 25 rows\n"
     ]
    }
   ],
   "source": [
    "def tcell_fraction(row):\n",
    "    return row[\"T-cell fraction\"]\n",
    "\n",
    "def peripheral_clonality_a(row):\n",
    "    return row['Clonality']\n",
    "\n",
    "cols, d = cohort.as_dataframe([snv_count,\n",
    "                               missense_snv_count,\n",
    "                               neoantigen_count,\n",
    "                               expressed_exonic_snv_count,\n",
    "                               expressed_missense_snv_count,\n",
    "                               expressed_neoantigen_count,\n",
    "                               exonic_snv_count,\n",
    "                               peripheral_clonality_a,\n",
    "                               tcell_fraction,\n",
    "                               ],\n",
    "                              rename_cols=True,\n",
    "                              return_cols=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:21.537104",
     "start_time": "2016-08-17T07:59:21.513702"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/plain": [
       "['snv_count',\n",
       " 'missense_snv_count',\n",
       " 'neoantigen_count',\n",
       " 'expressed_exonic_snv_count',\n",
       " 'expressed_missense_snv_count',\n",
       " 'expressed_neoantigen_count',\n",
       " 'exonic_snv_count',\n",
       " 'peripheral_clonality_a',\n",
       " 'tcell_fraction']"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "cols"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## construct/rescale variables"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:21.559557",
     "start_time": "2016-08-17T07:59:21.539192"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## add/modify count variables\n",
    "d['nonexonic_snv_count'] = d.snv_count - d.exonic_snv_count\n",
    "cols.append('nonexonic_snv_count')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:22.132709",
     "start_time": "2016-08-17T07:59:22.081884"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## create 'observed', log-transformed & centered versions of variables (not normalized by MB)\n",
    "for col in cols:\n",
    "    observed_col = 'observed_{}'.format(col)\n",
    "    log_col = 'log_{}'.format(col)\n",
    "    log_col_centered = 'log_{}_centered'.format(col)\n",
    "    log_col_rescaled = 'log_{}_rescaled'.format(col)\n",
    "    d[observed_col] = d[col]*d['mb']\n",
    "    d[log_col] = np.log1p(d[observed_col])\n",
    "    d[log_col_centered] = d[log_col] - np.mean(d[log_col])\n",
    "    d[log_col_rescaled] = d[log_col_centered]/np.std(d[log_col_centered])\n",
    "\n",
    "## save key vars in list for future use\n",
    "vars_centered = ['log_{}_centered'.format(col) for col in cols]\n",
    "vars_rescaled = ['log_{}_rescaled'.format(col) for col in cols]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:22.174374",
     "start_time": "2016-08-17T07:59:22.134744"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## construct new variables for key ratios / comparisons\n",
    "\n",
    "# what proportion of X are expressed?\n",
    "d['exonic_expression_ratio'] = d.expressed_exonic_snv_count / d.exonic_snv_count\n",
    "d['missense_expression_ratio'] = d.expressed_missense_snv_count / d.missense_snv_count\n",
    "d['neoantigen_expression_ratio'] = d.expressed_neoantigen_count / d.neoantigen_count\n",
    "\n",
    "# d['expressed_missense2snv_ratio'] = d.expressed_missense_snv_count / d.expressed_snv_count\n",
    "d['missense2exonic_snv_ratio'] = d.missense_snv_count / d.exonic_snv_count\n",
    "d['expressed_neoantigen2missense_ratio'] = d.expressed_neoantigen_count / d.expressed_missense_snv_count\n",
    "\n",
    "extra_cols = ['missense_expression_ratio','neoantigen_expression_ratio', 'exonic_expression_ratio', 'missense2exonic_snv_ratio','expressed_neoantigen2missense_ratio']\n",
    "## create recentered versions of ratios\n",
    "for col in extra_cols:\n",
    "    col_centered = '{}_centered'.format(col)\n",
    "    col_rescaled = '{}_rescaled'.format(col)\n",
    "    d[col_centered] = d[col] - np.mean(d[col])\n",
    "    d[col_rescaled] = d[col_centered]/np.std(d[col_centered])\n",
    "\n",
    "## append extra-cols to key var lists\n",
    "vars_centered.extend(['{}_centered'.format(col) for col in extra_cols])\n",
    "vars_rescaled.extend(['{}_rescaled'.format(col) for col in extra_cols])\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## center variables by mean within PD-L1 group"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:22.220909",
     "start_time": "2016-08-17T07:59:22.199735"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## identify list of variables to center\n",
    "metrics = list(cols)\n",
    "metrics.extend(extra_cols)\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "metrics2 = list(metrics)\n",
    "metrics2.extend(['pd_l1'])\n",
    "assert(not 'pd_l1' in metrics)\n",
    "log_metrics2 = ['log_{}'.format(var) for var in metrics2]\n",
    "metrics2.extend(log_metrics2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:22.242582",
     "start_time": "2016-08-17T07:59:22.222780"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "grp_metrics = [var for var in metrics2 if var in d.columns]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:22.275217",
     "start_time": "2016-08-17T07:59:22.244289"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# center variables by group\n",
    "bygrp = d.loc[:, grp_metrics]\n",
    "bygrp = bygrp.groupby('pd_l1').transform(lambda x: x - x.mean())\n",
    "bygrp['patient_id'] = d.patient_id"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:22.301832",
     "start_time": "2016-08-17T07:59:22.276693"
    },
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# merge recentered variables back into original dataframe\n",
    "df = pd.merge(d, bygrp, on = 'patient_id', suffixes = ['', '_centered_by_pd_l1'])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# prep data for survival analysis (long format)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T07:59:22.420738",
     "start_time": "2016-08-17T07:59:22.356618"
    },
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "## prep dflong_pfs which will be used for survival analysis using stan\n",
    "df_long_pfs = survivalstan.prep_data_long_surv(df = df,\n",
    "                                               event_col = 'is_progressed_or_deceased',\n",
    "                                               time_col = 'pfs')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# multivariate analysis using varying-coef model"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2016-08-17T08:27:09.524216",
     "start_time": "2016-08-17T08:27:09.501719"
    },
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "../utils/stan/logistic_model.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_custom_coefs_expressed_missense_and_neoant_mutations.stan\n",
      "../utils/stan/logistic_model_by_group.stan\n",
      "../utils/stan/pem_survival_model_unstructured_varcoef.stan\n",
      "../utils/stan/pem_survival_model_unstructured_varcoef_hsprior.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_tvc.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_alt.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_bspline_est_xi.stan\n",
      "../utils/stan/pem_survival_model_varying_coefs3.stan\n",
      "../utils/stan/pem_survival_model_randomwalk.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_custom_coefs_rate_only.stan\n",
      "../utils/stan/pem_survival_model_gamma.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_custom_coefs_missense_and_neoant_rates.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_tvc2.stan\n",
      "../utils/stan/pem_survival_model_varying_coefs2.stan\n",
      "../utils/stan/pem_survival_model_varying_coefs4.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_bspline_est_xi2.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_bspline.stan\n",
      "../utils/stan/pem_survival_model_randomwalk2.stan\n",
      "../utils/stan/pem_survival_model_randomwalk_bspline2.stan\n"
     ]
    }
   ],
   "source": [
    "models = survivalstan.utils.read_files('../utils/stan')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "survstan_pfs_varcoef = functools.partial(\n",
    "    survivalstan.fit_stan_survival_model,\n",
    "    df = df_long_pfs,\n",
    "    model_code = models['pem_survival_model_unstructured_varcoef.stan'],\n",
    "    timepoint_end_col = 'end_time',\n",
    "    event_col = 'end_failure',\n",
    "    sample_col = 'patient_id',\n",
    "    chains = 4,\n",
    "    iter = 10000,\n",
    "    grp_coef_type = 'vector-of-vectors',\n",
    "    seed = seed,\n",
    "    )\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## summarize varying-coef model for missense_snv_count by pd-l1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "NOT reusing model.\n"
     ]
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Ran in 161.237 sec.\n"
     ]
    },
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/home/tavi/miniconda2/lib/python2.7/site-packages/stanity/psis.py:228: FutureWarning: elementwise comparison failed; returning scalar instead, but in the future will perform elementwise comparison\n",
      "  elif sort == 'in-place':\n",
      "/home/tavi/miniconda2/lib/python2.7/site-packages/stanity/psis.py:246: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future\n",
      "  bs /= 3 * x[sort[np.floor(n/4 + 0.5) - 1]]\n"
     ]
    }
   ],
   "source": [
    "testfit = survstan_pfs_varcoef(formula = 'log_missense_snv_count_centered_by_pd_l1', \n",
    "                               group_col = 'pd_l1',\n",
    "                              )"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{{{hr_missense_by_pdl1_pfs}}}\n"
     ]
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAGHCAYAAACu4BXOAAAABHNCSVQICAgIfAhkiAAAAAlwSFlz\nAAALEgAACxIB0t1+/AAAIABJREFUeJzt3Xtc1HXe///nB0VQJyhELfFAhcUaKmbh5rKroXW1lZvo\narZGKm6paellpeVlWdd2aWltlq5rtaKpabqa6+Fa1xMewjS0SysULSVdCUXxgKCiiPP9w9/MDwTe\nzAwwg/G4327eos98Dq/Pez4zz/mc3h/LbrfbBQBAOfx8XQAAoGYjKAAARgQFAMCIoAAAGBEUAACj\nur4uwKSgoEBpaWlq3Lix6tSp4+tyAOC6UFRUpBMnTigqKkqBgYGVnl+NDoq0tDT179/f12UAwHXp\n008/1T333FPp+dTooGjcuLGkqyt78803+7ia619sbKwkKSUlxceVAKhOx44dU//+/Z3foZVVo4PC\ncbjp5ptvVvPmzX1czfXv0KFDvi4BgBdV1SF7TmYDAIwICgCAEUEBADAiKAAARgQFAMCIoKhFwsPD\nFR4e7usyAFxnCAoAgFGNvo8C1WfMmDHKycmpsvnl5+dLkmw2W5XN01WhoaGaPHmy15cL1BYERS2V\nk5OjE8ePq0EVze/C//df6/z5Kpqja7y7NKB2IihqsQaS+tStmk3g75cvS1U4P3eXC6D6cI4CAGDE\nHkUtQl9PADzBHgUAwIigAAAYERQAACOCAgBgRFAAAIwIilqEvp4AeIKgAAAYERQAACOCAgBgRFAA\nAIwICgCAEX091SL09QTAE+xRAACMCAoAgBFBAQAwIigAAEYEBQDAiKCooZKSkpSUlFSl86SvJ9RE\n1bGto2oRFDVUSkqKUlJSfF0GUO3Y1ms+ggIAYERQAACMCAoAgBFBAQAwoq+nWoS+ngB4gj0KAIAR\nQQEAMCIoAABGBAUAwIigAAAYERS1CH09AfAEQQEAMCIoAABGBAUAwMiloFi2bJkiIyN15MiREsNX\nrFihgQMHqlOnToqKilKXLl00evRopaamlhjvwIEDSkxMVIcOHdSpUye98sorys3Nrbq1AABUG5e7\n8LAsy/n3lStXNGrUKCUnJys+Pl4JCQkKDg7WsWPHtHr1ag0cOFCpqamy2Ww6fvy4EhISFBERoenT\npys3N1eTJ0/W0KFDtXDhwmpZKQBA1fGor6eZM2dq3bp1mjZtmrp3717itUcffVTbtm2Tv7+/JOlv\nf/ubioqK9Ne//lU2m02S1KRJEz355JNav359qelRfejrCYAn3D5HUVhYqDlz5qhr167lfsnfd999\nCggIkCRt3LhRXbp0cYaEJN1zzz1q1qyZNmzY4GHZAABvcTso0tLSdPbsWcXFxVU47sWLF5WZmanW\nrVuXei0iIkIHDhxwd/EAAC9zOyiOHj0qy7LUrFmzCsfNzc2V3W5XcHBwqdeCg4M5oQ0A1wEujwUA\nGLl9MvuWW26R3W5XVlZWheMGBQXJsqwy9xxyc3PL3NPAVfn5+SooKFBiYmK1zD8nJ0d1qmXO3nVJ\n0oWcnGprJ1S/nJwcBQYG+roMGLi9RxEVFaWgoCAlJydXOG5gYKDCwsLKPBdx4MABRUREuLt4VEJm\nZqYyMzN9XQaA64zbexT+/v4aNGiQPvjgA61du1YPPvhgqXG+/PJLdezYUQEBAYqLi9Py5cuVn5/v\nvPJp586dysrKUrdu3Sq/Bj9TNptNNptNSUlJVTZPR4eAa9euVWJios4dP15l8/aVepIahoZWaTvB\nu9gbrPk8uo9iyJAh2r9/v0aPHq2ePXvq/vvvV3BwsLKzs7VmzRqtX79eqampCggI0ODBg7Vy5UoN\nGzZMzzzzjPLy8vTOO+8oOjqaeygA4DrgUVD4+flp6tSpWrlypZYuXapx48bp3LlzCg0NVceOHTV/\n/nzn3kPTpk01d+5cvfXWW3r++edVr149devWTWPHjq3SFQEAVA+XgiI+Pl7x8fGlhvfo0UM9evSo\ncPrWrVtr1qxZ7lcHAPA5Lo8FABh5dOgJ1yf6egLgCfYoAABGBAUAwIigAAAYERQAACOCAgBgRFDU\nIuHh4c5uPADAVQQFAMCIoAAAGBEUAAAjggIAYERQAACM6OupFqGvJwCeYI8CAGBEUAAAjAgKAIAR\n5yhqqNjYWF+XAHgF23rNR1DUUImJib4uAfAKtvWaj0NPtQh9PQHwBEEBADAiKAAARgQFAMCIoAAA\nGBEUAAAjLo+tRejrCYAn2KMAABgRFAAAI4ICAGBEUAAAjAgKAIARQVGL0NcTAE8QFAAAI4ICAGBE\nUAAAjAgKAIARQQEAMKKvp1qEvp4AeIKgqMXOS/r75ctVNi9V4fzcWW5Dry4RqH0IiloqNDS0Sudn\nz8+XJDW02ap0vhVpqKpfFwAlERS11OTJk31dAoDrBCezAQBGBAUAwIigqEXo6wmAJwgKAIARQQEA\nMCIoAABGBAUAwIigAAAYccNdLUJfTwA8wR4FAMCIoAAAGBEUAAAjggIAYERQAACMCIpahL6eAHiC\noAAAGBEUAAAjggIAYERQAACMCAoAgBF9PXnZmDFjlJOTYxwnPz9fkmSz2SqcX2hoqCZPnuzSsunr\nCYAnCAovy8nJ0fHjJ1QnoH654xRdvCBJulhknpdjPACoTgSFD9QJqK8mHX9X7uvHv14hScZxio8H\nANWJcxQAACOCAgBgRFAAAIwIilqEvp4AeIKgAAAYERQAACOCAgBgRFAAAIwICgCAEXdm1yL09QTA\nE+xRAACMCAoAgBFBAQAwIigAAEYEBQDAiKCoRejrCYAnCAoAgBFBAQAwIigAAEYERSUkJSUpKSnJ\n12VUi5/zugFwD0FRCSkpKUpJSfF1GdXi57xuANxDX0+1CH09AfAEexQAACOCAgBgRFAAAIwICgCA\nEUEBADAiKGoR+noC4AmCAgBgRFAAAIwICgCAEUEBADAiKAAARvT1VIvQ1xMAT7BHAQAwcmmPYtmy\nZXrllVe0bt06tWjRwjl8xYoV+vzzz5Wenq5z586pUaNG6tixo/r166eYmBhJ0g8//KB58+YpLS1N\n33//vYqKipSenl49awMAqHIuH3qyLMv595UrVzRq1CglJycrPj5eCQkJCg4O1rFjx7R69WoNHDhQ\nqampstls2rNnj7744gtFRUUpICBAu3fvrpYVAQBUD4/OUcycOVPr1q3TtGnT1L179xKvPfroo9q2\nbZv8/f0lST179lTPnj0lSVOnTiUoAOA643ZQFBYWas6cOeratWupkHC47777Kl0YAKBmcPtkdlpa\nms6ePau4uLjqqAfViL6eAHjC7aA4evSoLMtSs2bNqqMeAEANw30UlZCfn6+CggIlJia6PE1OTo7s\nVtVclXzl8iXl5OS4vPw77rhDklwaPycnR4GBgZWqD8DPg9vfWLfccovsdruysrKqox4AQA3j9h5F\nVFSUgoKClJycrD59+lRHTdcNm80mm82mpKQkl6dJTEzUydz8Klm+X916ahTs+vId5yfWrl1b4bju\n7CUB+Hlze4/C399fgwYN0qZNm8r9wvnyyy918eLFShcHAPA9j85RDBkyRPv379fo0aPVs2dP3X//\n/QoODlZ2drbWrFmj9evXKzU1VQEBASooKNDmzZslSRkZGZKkNWvWSJLCwsIUFRVVRauCitDXEwBP\neBQUfn5+mjp1qlauXKmlS5dq3LhxOnfunEJDQ9WxY0fNnz9fNptNknTy5EmNHDmyxJ3do0aNknT1\nZrxJkyZVwWoAAKqLS0ERHx+v+Pj4UsN79OihHj16GKcNCwvTvn37PKsOAOBz9B4LADAiKAAARgQF\nAMCIoKhF6OsJgCcICgCAEUEBADAiKAAARgQFAMCIoAAAGPE8ilqEvp4AeII9CgCAEUEBADAiKAAA\nRgQFAMCIoAAAGBEUtQh9PQHwBEEBADDiPopKiI2N9XUJ1ebnvG4A3ENQVEJiYqKvS6g2P+d1A+Ae\nDj0BAIwICgCAEYeeahH6egLgCfYoAABGBAUAwIigAAAYERQAACOCAgBgRFDUIvT1BMATBAUAwIig\nAAAYERQAACOCAgBgRFAAAIzo66kWoa8nAJ5gjwIAYERQAACMCAoAgBFBAQAw4mS2DxRdvKDjX68w\nvi7JOM7/P56tKksDgFIICi8LDQ2tcJz8/Kv/tdkqCgGbS/NzcPTzxNVPANxBUHjZ5MmTfV0CALiF\ncxQAACOCAgBgRFAAAIwICgCAESezaxGudgLgCfYoAABGBAUAwIigAAAYERQAACOCAgBgRFDUIuHh\n4c7+ngDAVQQFAMCIoAAAGBEUAAAjggIAYERQAACM6OupFqGvJwCeYI8CAGBEUAAAjDj05EVjxoxR\nTk5Oua/n5+dLkmw2W7njhIaG8txtAF5FUHhRTk6Ojp84Lr/6ZTf7lQuXJUkFumR8HQC8iaDwMr/6\ndRXy21ZlvnZq9WFJqvB1APAmzlHUIvT1BMATBAUAwIigAAAYERQAACOCAgBgRFAAAIy4PLYWoa8n\nAJ5gjwIAYERQAACMCAoAgBFBAQAwIigAAEYERS1CX08APEFQAACMCAoAgBFBAQAwIigAAEYEBQDA\niL6eahH6egLgCfYoAABGBEUlJCUlKSkpyddlVMrPYR0AVC+CohJSUlKUkpLi6zIq5eewDgCqF0EB\nADAiKAAARgRFLUJfTwA8QVAAAIwICgCAEUEBADAiKAAARgQFAMCIvp5qEfp6AuAJ9igAAEYEBQDA\niKAAABgRFAAAI4ICAGDkUlAsW7ZMkZGROnLkSInhK1as0MCBA9WpUydFRUWpS5cuGj16tFJTU53j\nLFq0SIMHD1ZsbKyio6PVo0cPzZo1S4WFhVW7JqgQfT0B8ITLl8daluX8+8qVKxo1apSSk5MVHx+v\nhIQEBQcH69ixY1q9erUGDhyo1NRU2Ww2zZgxQ507d1afPn0UEhKir7/+Wu+//76+++47TZ06tVpW\nCgBQdTy6j2LmzJlat26dpk2bpu7du5d47dFHH9W2bdvk7+8vSfrHP/6hm266yfl6TEyMrly5ounT\npyszM1PNmzevRPkAgOrm9jmKwsJCzZkzR127di0VEg733XefAgICJKlESDi0bdtWkpSdne3u4gEA\nXuZ2UKSlpens2bOKi4vzeKGpqany8/PTrbfe6vE8AADe4XZQHD16VJZlqVmzZh4tcN++fZo3b556\n9+6tkJAQj+YBAPAer/b1dPz4cT377LNq1aqVXn75ZW8uulrk5+eroKBAiYmJLo2fk5OjK352j5d3\n5VKRcnJyXF7etRx7gcWnz8nJUWBgoMc1Afj5c3uP4pZbbpHdbldWVpZb0505c0aJiYny8/PTrFmz\n1KBBA3cXDQDwAbf3KKKiohQUFKTk5GT16dPHpWny8/OVmJio3NxcLViwQI0bN3a70JrIZrPJZrMp\nKSnJpfETExOVk3/K4+X51aujUFuIy8tztSYAMHF7j8Lf31+DBg3Spk2btHbt2jLH+fLLL3Xx4kVJ\nUkFBgZ555hllZWVp9uzZatGiReUqBgB4lUfnKIYMGaL9+/dr9OjR6tmzp+6//34FBwcrOztba9as\n0fr165WamqqAgACNGDFCu3fv1n/913/p3Llz+uabb5zzadGiBSe0AaCG8ygo/Pz8NHXqVK1cuVJL\nly7VuHHjdO7cOYWGhqpjx46aP3++bDabJCklJUWWZenNN98sNZ9JkyapZ8+elVsDAEC1ciko4uPj\nFR8fX2p4jx491KNHD+O0+/bt86wyVDlHP0886Q6AO+g9FgBgRFAAAIwICgCAEUEBADAiKAAARl7t\n6wm+xdVOADzBHgUAwIigAAAYERQAACOCAgBgRFAAAIwIilokPDzc2d8TALiKoAAAGBEUAAAjggIA\nYERQAACM6MKjEmJjY31dQqX9HNYBQPUiKCohMTHR1yW4pay+nq63dQDgfRx6AgAYERQAACOCAgBg\nRFAAAIwICgCAEUFRi9DXEwBPEBQAACOCAgBgRFAAAIwICgCAEUEBADCir6dapKy+ngCgIuxRAACM\nCAoAgBFBAQAwIigAAEYEBQDAiKCoRejrCYAnCAoAgBFBAQAw4oY7L7ty4bJOrT5c7muSzK/bqq00\nACgTQeFFoaGhxtfzlS9JstnKSQNbxfMAgKpGUHjR5MmTfV0CALiNoKhF6OsJgCc4mQ0AMCIoAABG\nBAUAwIigAAAYERQAACOCohahrycAniAoAABGNfo+iqKiIknSsWPHfFzJz0tmZqavSwBQjRzfmY7v\n0Mqq0UFx4sQJSVL//v19XMnPQ0BAgCSpW7duPq4EgDecOHFCrVq1qvR8LLvdbq+CeqpFQUGB0tLS\n1LhxY9WpU8fX5QDAdaGoqEgnTpxQVFSUAgMDKz2/Gh0UAADf42Q2AMCIoAAAGBEUAAAjggIAYOSz\ny2OPHTumiRMn6ssvv5Tdblfnzp01btw43XLLLRVOe+nSJb333ntauXKl8vLy9Itf/EIvvvii7rnn\nnhpTY2RkZKlhlmVp2bJlZb5WGdnZ2froo4+0Z88e7du3TwUFBUpOTlazZs0qnNZbbVnZOr3Vnv/6\n17+0cuVK7dmzR6dPn9Ytt9yiBx98UEOGDFHDhg2N03qzLStTpze3zZSUFH388cc6ePCgcnNzFRIS\nog4dOui5557T7bffbpzWm+1ZmTq92Z7XGjx4sLZu3aphw4Zp5MiRxnEr054+CYqCggI99dRTCggI\ncD717b333tOAAQO0YsWKCi/neuWVV/TFF19ozJgxat68uT799FMNHjxYixYtqrI3prI1SlLv3r31\n+OOPlxh26623Vkl9xR0+fFhr1qzRXXfdpXvuuUdbt251eVpvtGVV1Cl5pz1nz56tpk2b6oUXXtDN\nN9+s9PR0TZs2Tampqfrss8+M03qzLStTp+S9bTM3N1dRUVHq37+/QkJClJWVpY8++kiPP/64Vq5c\nafzR5c32rEydkvfas7hVq1Zp//79sizLpfEr1Z52H5gzZ469TZs29n//+9/OYUeOHLG3adPGPnv2\nbOO06enp9jvvvNO+bNky57DLly/b/+M//sM+bNiwGlGj3W6333nnnfapU6dWWT2uWrx4sT0yMtL+\n008/VTiut9qyLO7Uabd7rz1PnTpVatiyZcvskZGR9u3bt5c7nbfb0tM67XbfbZsOGRkZ9jvvvNP4\nOfLltungSp12u2/a88yZM/Zf/epX9v/93/91afmVbU+fnKPYuHGj2rdvrxYtWjiHNW/eXHfffbc2\nbNhgnHbDhg3y9/fXb3/7W+ewOnXq6JFHHlFKSooKCwt9XuP1wltteT256aabSg1r27at7Ha7srOz\ny53O223paZ01QXBwsCQZb6KtCdumK3X6yjvvvKM777xTDz/8sEvjV7Y9fRIUBw4cUOvWrUsNj4iI\n0MGDB43THjx4UM2bN3d2R1F82sLCQv373//2eY0OCxcuVNu2bRUdHa0BAwZo586dVVJbVfFWW1YV\nX7VnamqqLMsyHquuCW3pSp0O3m7LK1euqLCwUIcOHdKECRPUpEkTPfLII+WO76v2dLdOB2+2586d\nO7VixQq99tprLk9T2fb0yTmKM2fOONO6uODgYJ09e9Y4bW5ubpnT3njjjc55+7pGSXrsscfUtWtX\nNWnSRFlZWZo1a5YGDhyo2bNn6957762SGivLW21ZFXzVntnZ2Zo2bZo6d+6su+66q9zxfN2WrtYp\n+aYt+/Tpoz179kiSWrVqpTlz5igkJKTc8X3Vnu7WKXm3PQsLC/X6669r8ODBbvXhVNn2rNGdAl7P\n3n77beffHTt2VFxcnHr06KH3339f8+fP92Fl1ydftOf58+c1bNgw+fv7a+LEidWyjKrgbp2+aMsp\nU6YoPz9fmZmZmjVrlgYNGqSFCxe6dMWbN3lSpzfb8+OPP9bFixc1dOjQKp1vRXxy6Ck4OFi5ubml\nhufm5iooKMg4bVBQUJnTOhLRkZC+rLEsDRs2VJcuXfTdd99VRXlVwlttWR2quz0vXryoIUOG6Kef\nftKsWbPUtGlT4/i+akt36yyLN7bN2267Te3atdPDDz+sOXPm6Pz58/roo4/KHd9X7elunWWprvY8\nevSoPvzwQ40cOVIXL15UXl6e8+jGpUuXlJeXpytXrpQ5bWXb0ydBERERoQMHDpQafuDAgQqPr0ZE\nRCgzM1MXL14sNa2/v79atmzp8xqvF95qy+vN5cuX9dxzz2nv3r36+OOPFRERUeE0vmhLT+qsCW64\n4Qa1bNnSeFy8JmybrtTpTUeOHNGlS5f00ksv6d5779W9996rmJgYWZalWbNmKSYmRt9//32Z01a2\nPX0SFHFxcfrmm29KPEAnMzNTu3btqvBZCXFxcSosLNTq1audw4qKirR69WrFxsbK39/f5zWWJT8/\nX5s2bVK7du2qpL6q4K22rA7V1Z52u10vvPCCUlNTNWPGDJfn7+229LTOsnh728zJyVFGRobxy6km\nbJuu1FmW6mrPNm3aaO7cuZo7d67mzZvn/Ge32/XYY49p3rx55Z63qGx71nn99ddfr8qVccWdd96p\nf/7zn1qzZo2aNGmiH3/8URMmTFD9+vX15ptvOovOyspSp06dZFmW86RQ48aNlZGRoQULFujGG2/U\n2bNn9c477+i7777TO++8o9DQUJ/XmJSUpH/84x+6cOGCTp8+rdTUVI0fP15ZWVl6++23Xbqz211r\n1qzRwYMH9fXXXystLU3h4eHKysrS6dOnFRYW5tO2rGyd3mzP119/XcuXL9fTTz+tiIgIZWdnO/9Z\nliWbzVYj2tLTOr29bY4YMUKHDh3S2bNndeLECaWkpGjChAm6dOmSJk6cqBtvvLFGtKendXqzPevV\nq6ewsLBS/6ZPn664uDj16tVL/v7+1dKePjmZXb9+fX3yySeaOHGixo4d6+we45VXXlH9+vWd49nt\ndue/4t566y299957ev/995WXl6fIyEjNmjWrSu/WrEyNt956q9avX6+1a9cqLy9PNptNHTt21KRJ\nkxQVFVVlNRY3cuRI5x2almXpv//7vyVJ9957r+bOnevTtqxsnd5szy+++EKWZWnmzJmaOXNmideG\nDx+uESNG1Ii29LROb2+b0dHRWr16tebMmaPCwkLdfPPN6tSpk5555hnnCeKa0J6e1umLz/q1LMsq\ncXd2dbQnDy4CABjReywAwIigAAAYERQAACOCAgBgRFAAAIwICgCAEUEBADCi99hipk+frunTp0u6\neqfmiBEjSrzuuDHFsiylp6eXGu5gWZYaNmyo1q1bq2/fvoqPj3dp+YcOHdJbb72lb7/9VqdPn5bd\nbte4ceP01FNPVWa1XFJ83R38/PwUHBys9u3b649//GO1PKv452LatGn6y1/+IkmaN29ejelKvqr9\n9NNPJbqwuemmm7Rly5YSXUCkpaXp97//vfP/w8LCqu1hX45tNiwszOXPWVkSEhK0Y8eOEp/t1NRU\n52ev+PdBamqqUlNTJUm9evWqcT3gVgeCogymZ9CW99q1w8+dO6ddu3Zp165dOnfunJ588skKlztm\nzBh9++23znn5+Xl/h+/aOzzPnDmjTZs2acuWLfrggw/UvXt3r9d0PXDcHevq84uvd471PHPmjP71\nr3+pR48eztcWLFhQYpzq5AiKmJiYSgVFecp6X1NTUzV9+nRZlqVOnTrViqDg0JObyruR3TE8PT1d\nu3fv1vDhwyVd3dDmzZvn0rz37t0ry7J022236ZtvvtHevXurdG/C1cdHDh8+XOnp6dq5c6fzgfF2\nu71Ev/vluXTpUqVq9FRRUVG5XSx7w4gRI5Senq69e/de93sT7ryHdrtdCxcudP5/Xl6eVq9e7fxi\nra6OH4pvy9UVSDExMc731PF5rq0IimoQEBCggQMHSrr6QTl69Khx/GXLlikyMlJFRUWy2+06ePCg\n2rVrp8jISO3YsUOSdPr0aU2cOFEPPvig2rZtq7vvvlv9+vXT559/XmJeqampioyMVGRkpD744APN\nnDlTcXFxatOmjXbv3u3WejRs2FD/+Z//6VyPzMxMZ//1cXFxioyMVLdu3bRz507169dP7du314QJ\nE5zTf/7553riiSd09913q23btnrggQc0ceJEnT59usRyLl++rClTpig2NlYdOnTQH//4Rx0+fNi5\nHsUPdTjaKjIyUp999pneeustxcbGKioqSseOHZN09WlvEyZMULdu3RQVFaWYmBg9/fTTpR5Pefr0\nab3xxhvq3r27oqOj1bFjRz300EN64YUXdOjQIed4a9euVf/+/XXfffepbdu2io2N1ZNPPqnZs2c7\nx5k+fbqzLsd7Jl0NsDlz5qhXr17q0KGD2rVrp0ceeUQffPCBLly4UKIex/QJCQnavHmzevfurfbt\n2+uBBx7Q3/72N5ffN1fb3ZX30KRp06aqW7eudu3apR9++EHS1ffnwoULCgsLKzckqmJbXrVqlSIj\nI2VZlux2e4lxHT+utm3bpiFDhiguLk4dOnRQVFSUunbtqpdeesmlrsOLz9Ox5xIXF+fcm7Db7UpI\nSHCOk5KSotjYWEVGRurRRx8tMa8ffvjBOZ6r7VuTcOjJCxo1alThOMV/FRXvNE+62t1x3759lZWV\n5Rx2+fJl7d692/nP0ble8XksWLDA+bAST391mX6lW5alU6dOafDgwc5foY7lvPbaa1q8eHGJ5WZm\nZmru3LnasGGDFi9e7GyXCRMmaOnSpc5xt27dqoSEhAoPAU6dOrXU+mVkZOgPf/iDzpw54xyWl5en\nL774Qlu3btW7777rfMD82LFjtWXLlhLLOXz4sA4fPqzf/e53Cg8P17fffqtRo0aV+NI7efKkTp48\nqYKCAg0aNKhUXcXbbujQoc4O/BwyMjI0Y8YMbd68WZ9++qkCAwNLTJ+enl7iCWZHjhzRu+++q6ZN\nm5Y4xFMWd9rdsbxr30NXhYaGql27dlq7dq0WLlyo1157TZ999pksy9Ljjz+ud999t9Q0VbUtX3s4\nqKy/v/vuO23ZsqXEvLKzs7Vy5Upt27ZNq1atcukBSNduh+V9VgMCAtSvXz9Nnz5dBw8e1M6dO53n\n9VatWuV0d/clAAAKdElEQVQcr2/fvhUus6Zhj6IcxX8hOv65+mVbUFDg/LVpWVapXxfXio+PV3p6\nuux2u7Nr4OKHMaZOner8YPXq1UtfffWVli9f7jw2+ve//73MvYXc3FyNHz9eO3fu1MaNG3XHHXe4\n1Qb5+fmaOnWq8/9btGhR6oNVUFCgmJgYbdiwQbt27dLQoUP1f//3f84vq2bNmmn58uVKTU1Vr169\nJF3tmv3999+XdPUEviMkgoKCtGjRIn311Vfq0KFDmT1gOtjtdl24cEF//vOftWvXLq1Zs0YhISH6\nn//5H505c0ZBQUGaO3euvv32W61Zs0a33Xab7Ha7/vSnP+ny5cuSrj6k3rIsPfDAA9q5c6e+/vpr\nrVixQmPHjnU+Ke7rr792huWiRYuUlpamzZs3a+bMmRW+r6tWrXKGxC9+8QutX79eW7du1a9+9StJ\nVw81zp07t9R0586d09ChQ7Vjxw6NHz/eOXz58uXG5bnT7sVd+x4OGzbMuJzi+vXrJ0lasWKFNm3a\npIyMDAUGBjrPF1z7mamqbfn+++8v8zOTnp6uTz75RJIUGxur+fPna+vWrdqzZ4+++uorDRkyRNLV\nsF+xYoXL6+mQnJys4cOHO5c7b968Ep/VJ554QvXq1ZOkEofkHIfjIiMjK3ymeU1EUJSj+K8WV05S\nOnZFIyMjFR0drb/85S+qW7eu+vbtq5EjR1aqls2bNzv/Hjt2rIKCgnTHHXc4D29dO45D586d1b9/\nfzVs2FBNmzYt8+HqZXGE5D333KNFixZJunpifcyYMSXGc3yJT5o0Sc2aNVNgYKBatmxZopannnpK\nd9xxh2644Qa9/PLLznZ0/NLbvn27c9yePXuqXbt2CgoK0ujRoyWZLx547LHH9Nvf/laBgYFq0aKF\nLMvS9u3bZVmWzp49q4SEBLVt21YPPvigMjIyZLfbdfr0ae3du1eS1Lx5c9ntdu3atUszZszQmjVr\ndOnSJQ0YMMB5JVvz5s2dy/zwww/1ySefaO/evWrbtm2J9i9L8XZ49tlnFRYWppCQEL344otljuPQ\nqFEjPf/887LZbCVO0GZlZbm8vIra3aG899BVnTt3VqtWrXTu3DmNHTtWlmXpkUceKfdxwd7clps0\naaKVK1eqX79+io6OVkxMTIku2X/88UeX17M81/6QadSokR555BHZ7XatW7dOp0+f1nfffec81NWn\nT59KL9MXOPRUDkef/sW51G/7NVcNnT9/vtK1OI4tN2jQoMQHsPjVFidPniw1XZs2bTxaXvHd6aCg\nIEVHR2vw4MFlnqRt1KhRqYeenDp1qswab7jhBtlsNuXl5TnrLX7cvPhDXlx54Mu163fmzBkVFRUZ\ng92yLOcy33zzTb388sv68ccflZSU5PzQN2vWTDNmzFBkZKQeeOAB9e/fX0uWLFFycrKSk5Nlt9tV\np04d9evXT6+++mq59RVft+LtEBYW5vy7rPetZcuWzvobNGjgHH7tYyyv5U67F1fWe+iOvn37asqU\nKcrNzZVlWXriiSfKHddb27LdbteAAQN08ODBUodyHe9zQUGBW/N01YABA7Rs2TIVFhZqyZIlzvcl\nICBAv/vd76plmdWNPYoq4tgVTU9P16ZNmxQTE6OioiKtWrVKkydPrtS8Q0JCJEnnz59XXl6ec3jx\nk+RlnQcJCAjwaHmOq5727t2r7du3a+bMmeVeyVPWMhz1SiV/Befl5Sk/P1+WZTnrvemmm5yvZ2dn\nO/92nJg2KX5sX7r6gPg6depIklq1auU8FFH83969e9WlSxdJUrt27fTPf/5TGzZs0Mcff6wXX3xR\nDRo00NGjR/XOO+845/vqq69qx44dWrx4saZMmaIuXbqoqKhICxYs0DfffFNufeW1Q/G/y3rf6tb1\n7PebO+1enKfbiUOvXr1Ur149WZalu+66y3hoxVvb8v79+50hERERoY0bNyo9PV0zZsxwaz6eiIyM\nVExMjCRp8eLFzsNODz30kGw2W7UvvzoQFNWgadOmmjx5svNJeAsWLFBGRobH8+vatavz77fffltn\nz57V999/rzlz5pQ5jq85arHb7Zo3b56+//575eXl6a233nL+mnOM88tf/tI57ooVK7Rnzx7l5uY6\nT4S6c3llQECAfvnLX8put+vw4cOaMmWKTp06pcLCQmVkZGj27NkaMGCAc/z33ntPGzdulJ+fnzp1\n6qSHHnpIwcHBJa5U27Fjhz7++GNlZGQoPDxcDz74oNq3b++ch+lwUPH3ZObMmcrMzFROTk6JEKrK\n982ddq9KN910k4YPH65u3brp2WefdalGqWq25RtvvFF2u11ZWVk6e/asc7jjB4N09RGigYGB+umn\nn/Thhx+6PO/yFP9xs3///jK30QEDBjivFHT86Cl+E+L1hkNP1aRp06Z66qmn9OGHH6qoqEjvvfee\npk2bVuF0ZW10zz//vLZu3aqsrCwtWbJES5Yscb5mWZbzssaaokOHDnr88ce1ePFi/fTTTyV2ty3L\nUlhYmJ577jlJUnh4uH7/+99r6dKlOnnypHr37i3p6vHl4tO4aty4cerfv79yc3M1a9YszZo1q8Tr\nxQ/7rF69uswvDsuy9Otf/1rS1V+67777bplX8DRo0EAdO3Yst5aHH35YK1eu1JYtW5SWllbiZkXH\nr++EhIQS01R0n46JO+1eWdfW4zhJXNF4Vb0tR0dHa9OmTcrMzHT+ih8xYoSGDRum22+/XRkZGdqz\nZ4/zB0l4eHiZdbmjeH1vvvmm3nzzTUnSvn37nMPj4uLUsmVLHTlyRNLVPdzruWcD9iiuUdElmWUd\n/y5v+NNPP+28Smj9+vX69ttvK1x2WfMJDQ3V0qVLNWDAALVq1Ur16tVTw4YNFR0drUmTJpW6LtvT\nS2Hdnc50LuCNN97QpEmTFB0drYYNG8rf318tW7bUgAEDtGTJkhKHF15//XUNHjxYjRo1Uv369fWb\n3/xGH3zwgXMZ115pZVru7bffruXLl+uJJ55Qy5YtVa9ePQUFBal169bq06eP3njjDee4Tz75pO67\n7z41bdrU+auzdevWev755/XSSy9Jku666y717t1bERERCgoKUt26dRUSEqK4uDjNnTu3VKAVr8vP\nz09//etfNXbsWLVp00b169dXQECAIiIiNHz4cM2fP7/UpbHubF+VbfeK2rI87lzgce14Vb0tjx8/\nXl27dlVwcHCJ5dWpU0czZ87Ub37zG9lsNoWEhGjAgAEaP358he1c0fKjoqL06quvqmXLlvL395dl\nWaV6UbAsSwkJCc5D0tfjJbHF8cxs+NzBgwfl5+enW2+9VdLVk4yTJk3SokWLZFmWnnnmGeeNf8D1\n4s9//rM++ugjBQYGKjk5ucQ5pOsNh57gc9u3b9ef/vQnNWzYUEFBQcrJyVFhYaEsy9Ltt9+uwYMH\n+7pEwGVjxoxRamqqjh07Jsuy9Ic//OG6DgmJoEAN0KZNG/36179Wenq6cnJy5O/vr9atW6t79+4a\nOHBgiUtEgZru6NGjys7OVkhIiB5++GG98MILvi6p0jj0BAAw4mQ2AMCIoAAAGBEUAAAjggIAYERQ\nAACMCAoAgNH/A/W5Ez9OtvipAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<matplotlib.figure.Figure at 0x7fa98e3e2a10>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "paper.hyper_figure_label_printer('hr_missense_by_pdl1_pfs')\n",
    "#testfit['grp_coefs']['exp(beta)'] = np.exp(testfit['grp_coefs']['value'])\n",
    "sb.boxplot(data=testfit['grp_coefs'], x='exp(beta)', y='group', fliersize=0, whis=[2.5, 97.5])\n",
    "_ = plt.xlim([0, 4])\n",
    "_ = plt.vlines(1, -10, 10, linestyles='--')\n",
    "_ = plt.ylabel('')\n",
    "#_ = plt.title('Hazard associated with log(missense SNV count / MB) \\n by level of intratumoral PD-L1 expression')\n",
    "_ = plt.xlabel('HR for {}'.format(cohort.hazard_plot_name))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{{{hr_missense_pfs_among_IC0:n=4, HR=1.43, 95% CI (0.75, 2.98)}}}\n",
      "{{{hr_missense_pfs_among_IC1:n=11, HR=0.75, 95% CI (0.47, 1.14)}}}\n",
      "{{{hr_missense_pfs_among_IC2:n=10, HR=0.73, 95% CI (0.48, 1.06)}}}\n"
     ]
    }
   ],
   "source": [
    "testfit['grp_coefs']['group_label'] = testfit['grp_coefs'].group.apply(\n",
    "    lambda x: 'hr_missense_pfs_among_{}'.format(x)\n",
    ")\n",
    "for (name, ic_class), group in testfit['grp_coefs'].groupby(['group_label','group']):\n",
    "    group_n = len(df.loc[df.pd_l1==ic_class,'pd_l1'])\n",
    "    paper.hyper_label_printer(formatter=paper.hr_posterior_formatter, label=name, n=group_n, series=group['exp(beta)'])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "{{{pval_missense_pfs_interaction:0.046}}}\n"
     ]
    }
   ],
   "source": [
    "## bayesian p-value for interaction\n",
    "comparison = pd.pivot_table(testfit['grp_coefs'],\n",
    "                  index = ['iter', 'model_cohort', 'variable'],\n",
    "                  values = 'value', columns = 'group')\n",
    "print('{{{%s:%s}}}' % ('pval_missense_pfs_interaction', paper.float_str(1-(comparison.eval('IC2 < IC0').mean()))))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 2",
   "language": "python",
   "name": "python2"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.11"
  },
  "name": "TCR Clonality.ipynb",
  "nav_menu": {},
  "toc": {
   "navigate_menu": true,
   "number_sections": true,
   "sideBar": true,
   "threshold": 6,
   "toc_cell": false,
   "toc_section_display": "block",
   "toc_window_display": true
  },
  "toc_position": {
   "height": "877px",
   "left": "0px",
   "right": "2073px",
   "top": "106px",
   "width": "538px"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}
